module ice_probability

    use ice_kinds_mod 
    use ice_communicate, only: my_task, master_task, get_num_procs
    use ice_constants, only: pi, puny, c1
    use ice_fileunits, only: nu_timing
    use ice_domain_size, only: nx_global, ny_global, max_blocks
    use ice_domain, only: nblocks, blocks_ice, distrb_info, distribution_wght_file
    use ice_blocks, only: nx_block, ny_block,  nblocks_tot, block, get_block
    use ice_read_write, only: ice_open_nc, ice_read_global_nc, ice_close_nc   
    use ice_broadcast, only: broadcast_array
    use ice_global_reductions, only: sum_vector_dbl
    use ice_timers, only: 
    use ice_work, only: work_gr
    use ice_probability_tools

    implicit none 

   integer (int_kind), public, parameter ::  &    ! types of blocks:
         lndType     = 0,                    &    !     Land
         icefreeType = 1,                    &    !     ice free (ocean only)
         iceType     = 2                          !     sea ice 

!    include "netcdf.inc"
   integer (int_kind), public              :: dynCnt          ! number of calls to step_dynamic
   real (dbl_kind), allocatable :: lnumIceCells(:)  ! number of active ice cells

   public :: ReadProbabilityFile

   public :: CalcWorkPerBlock

!   public :: WriteProbabilityStats

   public :: init_numIceCells,   &
             accum_numIceCells,  &
             accum_numIceCells2, &
             print_numIceCells

   public :: set_numIceCells,write_numIceCells 
   

contains 
   
   subroutine init_numIceCells()

       dynCnt = 0

!       print *,'init_numIceCells: nblocks_tot is: ',nblocks_tot
       allocate(lnumIceCells(nblocks_tot))
       lnumIceCells = 0.0

   end subroutine init_numIceCells

   subroutine accum_numIceCells(iblk,icells)

     integer (int_kind) :: iblk,icells

     lnumIceCells(iblk) = lnumIceCells(iblk) + real(icells,kind=dbl_kind)

   end subroutine accum_numIceCells

   subroutine accum_numIceCells2(aice)

     real (dbl_kind) :: aice(nx_block,ny_block,max_blocks)
 
     integer (int_kind) :: igblk,iblk

     real (dbl_kind) :: tmp
     type (block) :: this_block
     integer (int_kind) :: i,j,ihi,ilo,jhi,jlo
  
  
     do iblk = 1,nblocks     
        igblk = blocks_ice(iblk)
        this_block = get_block(igblk,iblk)
        ilo = this_block%ilo
        ihi = this_block%ihi
        jlo = this_block%jlo
        jhi = this_block%jhi
        tmp = 0.0
        do j = jlo,jhi
        do i = ilo,ihi
           if(aice(i,j,iblk) > puny) then 
              tmp = tmp + c1
           endif
        enddo
        enddo 
        lnumIceCells(igblk) = lnumIceCells(igblk) + tmp
     enddo

   end subroutine accum_numIceCells2



   subroutine set_numIceCells(iblk,ncells)
      integer (int_kind) :: iblk
      real (dbl_kind)  :: ncells
    
     lnumIceCells(iblk) = ncells 
  
   end subroutine set_numIceCells

   subroutine write_numIceCells
    real (dbl_kind), allocatable :: gnumIceCells(:)
    integer :: n

    allocate(gnumIceCells(nblocks_tot))
    
    call sum_vector_dbl(lnumIceCells,gnumIceCells,distrb_info)

    gnumIceCells=gnumIceCells/real(dynCnt,kind=dbl_kind)
    if(my_task == master_task) then
       open(nu_timing,file='numCells2.bin',recl=8*nblocks_tot, &
          form = 'unformatted', access = 'direct', status = 'unknown')
        write(nu_timing,rec=1) gnumIceCells
        close(nu_timing)
    endif

 

   end subroutine write_numIceCells

   subroutine print_numIceCells

      real (dbl_kind), allocatable :: gnumIceCells(:)
      real (dbl_kind), allocatable :: gnumIceCells2(:)
      integer :: ii,n

      allocate(gnumIceCells(nblocks_tot))
      allocate(gnumIceCells2(nblocks_tot))

      call sum_vector_dbl(lnumIceCells, gnumIceCells, distrb_info)

      gnumIceCells = gnumIceCells/real(dynCnt,kind=dbl_kind)
      if(my_task == master_task) then 
         ! compress out land blocks 
         ii =0
         do n=1,nblocks_tot
            if(nocn(n) > 0) then 
               ii = ii+1
               gnumIceCells2(ii) = gnumIceCells(n)
            endif
         enddo
          open(nu_timing,file='numCells.bin',recl=8*ii, &
            form = 'unformatted', access = 'direct', status = 'unknown')
          write(nu_timing,rec=1) gnumIceCells2(1:ii)
          close(nu_timing)
          print *,'numCells: ',gnumIceCells2(1:ii)
      endif

      deallocate(gnumIceCells,gnumIceCells2)
    

   end subroutine print_numIceCells

   subroutine ReadProbabilityFile(distribution_wght_file,Prob)
   
      character(char_len_long), intent(in) :: distribution_wght_file
      real(real_kind), intent(inout)           :: Prob(:,:)
      
      type(block) :: this_block
      integer(int_kind) :: ilo,ihi
      integer(int_kind) :: jlo,jhi

      integer(int_kind) :: fid_prob
      integer(int_kind) :: amode,ncid,iostat,varid
      integer(int_kind) :: i,j,n,ncnt,ierr,ig,jg
      real(dbl_kind) :: val,sum,avg

      character (char_len) :: &
         fieldname       ! field name in netCDF file


      call ice_open_nc(distribution_wght_file,fid_prob)


      fieldname = 'ice_present'
      call ice_read_global_nc(fid_prob,1,fieldname,Prob,.true.)

      if(my_task == master_task) then 
        call ice_close_nc(fid_prob)
        print *,'MAXVAL(Prob): ',MAXVAL(Prob) 
        print *,'MINVAL(Prob): ',MINVAL(Prob) 
        print *,'COUNT(Prob > 0.5): ',COUNT(Prob > 0.5)
      endif

      ! normalize probability [0,1] 
!      prob=prob/100.0
      
   end subroutine ReadProbabilityFile

!***********************************************************************

!***********************************************************************
 subroutine CalcWorkPerBlock(distribution_wght, KMTG,ULATG,work_per_block, prob_per_block,blockType,bStats) 

   character (char_len), intent(in) :: distribution_wght

   real (dbl_kind), dimension(nx_global,ny_global), intent(in) :: &
      KMTG           ,&! global topography
      ULATG            ! global latitude field (radians)

   integer (int_kind), intent(inout)   :: &
      work_per_block(nblocks_tot)         ! number of work units per block

   real (dbl_kind), intent(inout) :: &
      prob_per_block(nblocks_tot)       ! probability of sea-ice in block

   integer (int_kind), intent(inout) ::  &
      blockType(nblocks_tot)            ! Type of block: 
                                        ! one of the following
                                        !    lnd,ice,icefree

   real (dbl_kind), intent(inout) :: bStats(:,:) 


   type (block) :: &
      this_block           ! block information for current block

   integer (int_kind), parameter :: &
      max_work_unit=10     ! quantize the work into values from 1,max

   integer(int_kind) :: jg,j,n,ig,i,work_unit
   integer(int_kind) :: maxkmt,mkmt,ilo,jlo,ihi,jhi
   integer(int_kind) :: activePts

   integer(int_kind) :: bid(max_blocks)
   integer(int_kind) :: bsize_x,bsize_y

   integer(int_kind), allocatable :: iLocation(:)

   real(dbl_kind), allocatable, dimension(:)   :: prob,work

   real(dbl_kind) :: lat
   real(dbl_kind), parameter :: w0 = 1., &  ! Constants for work blocks:
                          w1 = 10.   ! w0 for all blocks
                                      ! w1 for sea-ice blocks
   real(dbl_kind) :: maxlat,maxalat
   real(dbl_kind) :: shlatT,nhlatT
   real(real_kind) :: tmp

   integer(int_kind) :: ActiveBlocks,totWork

   logical, parameter :: Debug = .FALSE.

   integer(int_kind) :: numLnd,numIcefree,numIce

   integer(int_kind) :: ierr
   
!----------------------------------------------------------------------
!
!  estimate the amount of work per processor using the topography
!  and latitude
!
!----------------------------------------------------------------------



   allocate(work_gr(nx_global,ny_global)) 
   allocate(prob(nblocks_tot),work(nblocks_tot))
   allocate(nocn(nblocks_tot))
   allocate(nice005(nblocks_tot),nice010(nblocks_tot),nice050(nblocks_tot), &
	    nice100(nblocks_tot),nice250(nblocks_tot),nice500(nblocks_tot))

   nocn=0
   nice005=0
   nice010=0
   nice050=0
   nice100=0
   nice250=0
   nice500=0

   ActiveBlocks = 0
   select case (distribution_wght)
     case('file')
        call ReadProbabilityFile(distribution_wght_file,work_gr)
     case('erfc')
        work_gr = ErfcProbability(ULATG)
     case default
        work_gr = ErfcProbability(ULATG)
   end select
   !--------------------------------------------------
   ! It would be nice if we can not find the 
   ! the distribution_wgt_file,
   ! fall back to erfc weight function.
   ! However currently the error trapping does not prevent 
   ! such a fault recovery 
   !-------------------------------------------------------
   !if(ierr < 0) then 
   !    print *,'CalcWorkPerBlock: Could not open file: ',trim(distribution_wght_file)
   !    print *,'CalcWorkPerBlock:  Using ERFC probability function instead'
   !    work_gr = ErfcProbability(ULATG)
   !endif
       
     
!   print *,'IAM: ',my_task,'prod: ',prob

if(my_task == master_task) then 
   work_per_block=0
   do n=1,nblocks_tot

      this_block=get_block(n,n)
      ilo = this_block%ilo;ihi = this_block%ihi
      jlo = this_block%jlo;jhi = this_block%jhi
      !----------------------------------------------------
      ! calculate the probability of sea-ice in this block 
      !----------------------------------------------------
      tmp = 0.0
      do j=jlo,jhi
         jg=this_block%j_glob(j)
         if(jg>0) then  
            do i=ilo,ihi
              ig = this_block%i_glob(i)
              if(ig>0) then 
                 if(KMTG(ig,jg)>puny) then 
                    nocn(n) = nocn(n) + 1
                    if(work_gr(ig,jg) > 0.005) nice005(n) = nice005(n) + 1  
                    if(work_gr(ig,jg) > 0.010) nice010(n) = nice010(n) + 1  
                    if(work_gr(ig,jg) > 0.050) nice050(n) = nice050(n) + 1  
                    if(work_gr(ig,jg) > 0.100) nice100(n) = nice100(n) + 1  
                    if(work_gr(ig,jg) > 0.250) nice250(n) = nice250(n) + 1  
                    if(work_gr(ig,jg) > 0.500) nice500(n) = nice500(n) + 1  
                    tmp = tmp + work_gr(ig,jg)
                 endif
              endif
           enddo
         endif
      enddo
      if(nocn(n) > 0) then 
!         print *,'n:',n,' tmp:',tmp,' nocn(n): ',nocn(n)
         prob_per_block(n) = real(tmp,kind=dbl_kind)/real(nocn(n),kind=dbl_kind)
!         print *,'prob_per_block(n):' ,prob_per_block(n)
!         print *,' ihi,ilo,(ihi-ilo+1): ',ihi,ilo,(ihi-ilo+1)
!         print *,' jhi,jlo,(jhi-jlo+1): ',jhi,jlo,(jhi-jlo+1)
!         print *,' w0,w1: ',w0,w1
         work_per_block(n) =  &
           CEILING(w0 + (real(nocn(n),kind=dbl_kind)/real((ihi-ilo+1)*(jhi-jlo+1),kind=dbl_kind)) & 
                            *prob_per_block(n)*w1,kind=int_kind)
      else
         prob_per_block(n) = 0.0d0
      endif

      !--------------------------------------------------
      ! set the type of block (used latter for partition) 
      !--------------------------------------------------
      if(prob_per_block(n) >0.005 .and. nocn(n) > 0) then
        blockType(n) = iceType
      elseif (nocn(n) > 0) then 
        blockType(n) = icefreeType
      elseif (nocn(n) == 0) then
        blockType(n) = lndType
      endif

   enddo
 
   numIceFree = COUNT(blockType .eq. iceFreeType) 
   numLnd     = COUNT(blockType .eq. lndType)
   numIce     = COUNT(blockTYpe .eq. iceType)
    
!   print *,'Total blocks:',nblocks_tot,' land blocks: ',numLnd,' Ice blocks: ', &
!		numIce,' IceFree blocks: ',numIceFree
   write(*,23) nblocks_tot,numIce,numIceFree,numLnd 

23   format('CalcWorkPerBlock: Total blocks: ',i5,' Ice blocks: ',i5,' IceFree blocks: ',i5,' Land blocks: ',i5)

endif
   !-----------------------------------------------------------------
   ! distribute blocks among processors
   !-----------------------------------------------------------------
   call broadcast_array(work_per_block,master_task)
   call broadcast_array(prob_per_block,master_task)
   call broadcast_array(blockType,master_task)
   call broadcast_array(nocn,master_task)
   call broadcast_array(nice005,master_task)
   call broadcast_array(nice010,master_task)
   call broadcast_array(nice050,master_task)
   call broadcast_array(nice100,master_task)
   call broadcast_array(nice250,master_task)
   call broadcast_array(nice500,master_task)

   allocate(iLocation(nblocks_tot))
   do i=1,nblocks_tot
       iLocation(i) = i
   enddo

   call BuildProbabilityStats2(iLocation,bStats)

   deallocate(work_gr)
   deallocate(iLocation) 
!DBG   print *,'CalcWorkPerBlock: at the end of the subroutine' 
      
 end subroutine CalcWorkPerBlock

 function ErfcProbability(lat) result(prob)

      real (kind=dbl_kind), intent(in), &
         dimension (nx_global,ny_global) :: lat

      real (kind=real_kind),  &
	 dimension (nx_global,ny_global) :: prob


      real (kind=dbl_kind) :: ltmp,thetai,sigmai,arg
      integer :: i,j

      
      do j=1,ny_global
      do i=1,nx_global
         ltmp = lat(i,j)         
         if(ltmp > 0.0d0) then 
                   !---------------------
           ! northern latitude
           !---------------------
!JMD           thetai = (70.0_dbl_kind/180.0_dbl_kind)*pi
!JMD  For startup lots of random ice at low lattitudes
           thetai = (55.0_dbl_kind/180.0_dbl_kind)*pi
           sigmai = (5.0_dbl_kind/180.0_dbl_kind)*pi
           arg = (thetai-ltmp)/sigmai
        else 
           !---------------------
           ! southern hemisphere
           !---------------------
!JMD           thetai = (60.0_dbl_kind/180.0_dbl_kind)*pi
!JMD  For startup lots of random ice at low lattitudes
           thetai = (60.0_dbl_kind/180.0_dbl_kind)*pi
           sigmai = (5.0_dbl_kind/180.0_dbl_kind)*pi
           arg = (thetai+ltmp)/sigmai
	 endif
         prob(i,j) = real(erfc06(arg),kind=real_kind)/2.0_real_kind
      enddo
      enddo

 end function ErfcProbability 

 function erfc06(x)  result(erfc)

      implicit none
      real(dbl_kind), intent(in)     :: x
      real(dbl_kind)                 :: erfc
      real(dbl_kind)                 :: t,u
      real(dbl_kind), parameter      :: pa  =  3.97886080735226000e+00_dbl_kind
      real(dbl_kind), parameter      :: p00 =  2.75374741597376782e-01_dbl_kind
      real(dbl_kind), parameter      :: p01 =  4.90165080585318424e-01_dbl_kind
      real(dbl_kind), parameter      :: p02 =  7.74368199119538609e-01_dbl_kind
      real(dbl_kind), parameter      :: p03 =  1.07925515155856677e+00_dbl_kind
      real(dbl_kind), parameter      :: p04 =  1.31314653831023098e+00_dbl_kind
      real(dbl_kind), parameter      :: p05 =  1.37040217682338167e+00_dbl_kind
      real(dbl_kind), parameter      :: p06 =  1.18902982909273333e+00_dbl_kind
      real(dbl_kind), parameter      :: p07 =  8.05276408752910567e-01_dbl_kind
      real(dbl_kind), parameter      :: p08 =  3.57524274449531043e-01_dbl_kind
      real(dbl_kind), parameter      :: p09 =  1.66207924969367356e-02_dbl_kind
      real(dbl_kind), parameter      :: p10 = -1.19463959964325415e-01_dbl_kind
      real(dbl_kind), parameter      :: p11 = -8.38864557023001992e-02_dbl_kind
      real(dbl_kind), parameter      :: p12 =  2.49367200053503304e-03_dbl_kind
      real(dbl_kind), parameter      :: p13 =  3.90976845588484035e-02_dbl_kind
      real(dbl_kind), parameter      :: p14 =  1.61315329733252248e-02_dbl_kind
      real(dbl_kind), parameter      :: p15 = -1.33823644533460069e-02_dbl_kind
      real(dbl_kind), parameter      :: p16 = -1.27223813782122755e-02_dbl_kind
      real(dbl_kind), parameter      :: p17 =  3.83335126264887303e-03_dbl_kind
      real(dbl_kind), parameter      :: p18 =  7.73672528313526668e-03_dbl_kind
      real(dbl_kind), parameter      :: p19 = -8.70779635317295828e-04_dbl_kind
      real(dbl_kind), parameter      :: p20 = -3.96385097360513500e-03_dbl_kind
      real(dbl_kind), parameter      :: p21 =  1.19314022838340944e-04_dbl_kind
      real(dbl_kind), parameter      :: p22 =  1.27109764952614092e-03_dbl_kind




      t = pa/(pa + abs(x))
      u = t - 0.5_dbl_kind
      erfc = ((((((((((((((((((((((p22 * u + p21) * u + p20)   &
                                       * u + p19) * u + p18)   &
                                       * u + p17) * u + p16)   &
                                       * u + p15) * u + p14)   &
                                       * u + p13) * u + p12)   &
                                       * u + p11) * u + p10)   &
                                       * u + p09) * u + p08)   &
                                       * u + p07) * u + p06)   &
                                       * u + p05) * u + p04)   &
                                       * u + p03) * u + p02)   &
                                       * u + p01) * u + p00) * t * exp(-x**2)
      if (x < 0.0_dbl_kind) erfc = 2.0_dbl_kind - erfc

      end function erfc06

!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||

end module ice_probability
